
##########################################################
## This code is to pack the MODIS granule AOD date into ##
## daily datafile.                                      ##
##########################################################
rm(list=ls())
options(scipen=100,digits=10)

## Note: rgdal and gdalUtils packages are only available in the prevous R versions (R 4.2.0 and earlier)
library(rgdal)
library(gdalUtils)
library(raster)
library(dplyr)
library(lubridate)


setwd("raw data/AOD") 
# This is the raw satellite data downloaded from NASA (https://search.earthdata.nasa.gov/)  We use MOD04_3K from Terra. 
out.files <- list.files(getwd(), pattern="hdf$", full.names=FALSE) #create a list with names of the .hdf files (they should be stored in your workspace)
file.name=unique(regmatches(out.files, regexpr("A.......", out.files)))
date.name=unique(regmatches(file.name, regexpr("[0-9].*[0-9]", file.name)))

# Due to the large data size, we process this data by granule.
for(i in 1:length(date.name)){
  setwd("raw data/AOD")
  files=list.files(getwd(), pattern=file.name[i], full.names=FALSE) 
  for(n in 1:length(files)){
    # Read HDF file
    a=get_subdatasets(files[n])
    # Translate HDF file layer 16 to readable Raster dataframe. 
    # Layer 16: Corrected land AOD
    gdal_translate(a[16], dst_dataset = "AOD.tif")
    rast_AOD <- raster("AOD.tif")
    # Transfer raster dataframe to spatial dataframe
    AOD=rasterToPoints(rast_AOD, spatial=TRUE)
    # Transfer spatial dataframe to viewable dataframe to extract AOD value
    AOD=as.data.frame(AOD)
    
    ## No.52 layer is Longitude 
    gdal_translate(a[52], dst_dataset = "long.tif")
    rast_long <- raster("long.tif")
    Long=rasterToPoints(rast_long, spatial=TRUE)
    # Transfer spatial dataframe to viewable dataframe to extract longitude value
    Long=as.data.frame(Long)
    
    ## No.53 layer is Latitute 
    gdal_translate(a[53], dst_dataset = "lat.tif")
    rast_lat <- raster("lat.tif")
    Lat=rasterToPoints(rast_lat, spatial=TRUE)
    # Transfer spatial dataframe to viewable dataframe to extract latitude value
    Lat=as.data.frame(Lat)
    
    # Create dataframe with AOD, Long and Lat.
    AOD_Geo=left_join(AOD, Long,by=c("x","y") )
    AOD_Geo=left_join(AOD_Geo,Lat,by=c("x","y"))
    AOD_Geo=AOD_Geo[c(1,4,5)]
    # Store AOD data for each granule
    assign(paste("AOD_Geo_", n, sep = ""), AOD_Geo)
  }
  
  AOD_Geo_full=data.frame(matrix(ncol = 3, nrow = 2))
  x <- c("long", "lat", "AOD")
  colnames(AOD_Geo_full) <- x
  
  for(nn in 1:length(files)){
    # Get the name of granule data
    nam = paste('AOD_Geo_',nn,sep='')
    # Get the granule data
    AOD_temp=get(nam)
    # Put the granule data into dailty data file
    AOD_Geo_full=full_join(AOD_Geo_full,AOD_temp,by=c("long","lat"))
    # Give same name to same columns for next step of merging two columns
    names(AOD_Geo_full)[3]="AOD"
    names(AOD_Geo_full)[4]="AOD"
    
    # Merge two AOD columns into single column
    
    AOD_Geo_full <- as.data.frame( sapply(unique(names(AOD_Geo_full)), # for each unique column name
                                          function(col) rowMeans(AOD_Geo_full[names(AOD_Geo_full) == col],na.rm = TRUE) # calculate row means
    )
    )
    
    
  }
  AOD_Geo_full=AOD_Geo_full[-(1:2),]
  rownames(AOD_Geo_full) <- NULL
  
  # Set the date variable
  day=substr(date.name[i] , 5,7)
  year=substr(date.name[i] , 1,4)
  origin <- as.Date(paste0(year, "-01-01"),tz = "UTC") - days(1)
  date=as.Date(as.numeric(day), origin = origin, tz = "UTC")
  if(nrow(AOD_Geo_full)!=0) {
    AOD_Geo_full$Time=as.Date(date)
    # Store the processed daily AOD datafile (data size too large, did not included in the replication package)
    setwd("processed data/AOD_processed")
    write.csv(AOD_Geo_full,paste('AOD_MOD04_3K_',date,sep='',".csv"))
  }
  
  
  
  
}
